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ABSTRACT 



Using linear systems theory as a framework, the solution for the acoustic field 
present in a range-independent acoustic channel excited by a complex -weighted, 
planar array of point sources with an arbitrary input electrical signal is derived. 
The ocean medium is characterized by a transfer function, obtainable as the solution 
to the Helmholtz wave equation. The transfer function for an isospeed, three-layer 
waveguide is derived. The unbounded homogeneous medium equations are derived 
as a special case of the waveguide problem. The problem of interference due to the 
presence of a pressure -release surface is also derived as a special case. The linear 
systems approach lends itself to a modular computer implementation, in which 
different ocean medium models are represented by subroutines implementing their 
transfer functions. The equations for a range-independent medium are implemented 
as a group of subprograms. Results are presented for the special cases of a 
homogeneous medium and the surface reflection problem, which can be checked 
against known, easily interpreted analytical solutions. Finally, an example of 
waveform prediction for the isospeed, three-layer waveguide is presented. 
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I. INTRODUCTION 



Model-based or matched-field signal processing refers to the use of the 
knowledge of the physics of a problem for the construction of mathematical models 
and its application to suitable signal processing algorithms. Its application to 
underwater acoustics has been the subject of many papers in the signal processing 
literature, as in Ziomek and Blount [Ref. 1], Baggeroer, Kuperman and 
Schmidt [Ref. 2] and references therein. As the interest in this field grows and new 
algorithms are developed, so will the need for simulation tools to verify the behavior 
of those algorithms. Such tools should be able to generate signals with spatial and 
temporal structures like those found in real environments. 

The waveform prediction problem has been studied by Officer [Ref. 3:pp. 
101-110.130] for the specific cases of pulse reflection from a boundary and 
transmission in shallow water. DiNapoli and Deavenport [Ref.4:pp. 131-134] 
describe briefly a method employing the Fast -Field -Program (FFP) suitable to 
range -independent problems. Our purpose is to derive a generic solution for the 
range -independent, time -invariant , deterministic acoustic channel excited by an 
array of point sources transmitting an arbitrary signal. In our derivation we will 
use a linear systems approach to the acoustic problem, as described in 
Ziomek [Ref. 5], whose notation we follow closely. In particular, we will apply our 
results to the isospeed Pekeris waveguide. As we will show, the linear systems 
approach lends itself to a modular computer implementation, where different ocean 
media are represented by subroutines or functions implementing their transfer 
function. 
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In Section II we present the basic input -output relationships for a 
space-variant, time-independent linear system, which are the basis of our 
derivations. Next, we show how those equations can be used to represent an 
acoustic channel. The particular case of a range -independent ocean is then studied, 
when we derive the output waveform equation and describe the process of derivation 
of the medium transfer function. Finally, we apply the above results to the isospeed 
Pekeris waveguide. In Section III we discuss computer implementation issues and 
present some results of a computer implementation of the Pekeris waveguide 
equations. 
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II. ANALYSIS 



A. THE ACOUSTIC CHANNEL AS A LINEAR SYSTEM 

In this section we will present some results from the theory of space -time linear 
systems [Ref. 5], establish the notation to be used, and derive some intermediate 
results to be used in subsequent sections. 

1. Space - Vari cm t, Time-Invariant Linear Systems 
a. Impulse Response 

Linear, time-invariant, space-variant filters, as shown in Fig. 1, are 
represented by linear partial differential equations whose coefficients are time 
independent. In linear systems terminology, such systems are characterized by their 
time-invariant, space -dependent impulse response h(r, r 0 ; r). The output ^(i.r) is 
given by [Ref. 5;p. 8] 



«'-J7 



g(t-T, r-r 0 ) A(r,r 0 ;r) dr dr c 



( 1 ) 





A(r,r 0 ;r) 


t 


£(/.<-) 




*(/,») 



Fig. 1. A Linear, Time-Invariant, Space-Variant System 
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A physical interpretation of h, r and r 0 can be obtained by letting the 
input be an impulse applied at time t = to and at position r = r s , that is, 
g(£,r) = 6(t-t 0) T-i s ) = 6(t-to ) 6(r-r s ). Substituting this expression for g(t, r) into 
Eq. (1) and using the sifting property of the impulse function we obtain, for the 
output, 



p(f,r) = h( t-t 0 , r-r s ; r ), 

XX XX 

T r c 



(2) 



from which we can define A(r,r 0 ;r) as the response of the system at time t and 
position r due to the application of an impulse r seconds ago, that is, at the instant 
to = t-r, when the point source is positioned at r s = r-r 0 . Figure 2 illustrates the 
geometry of this problem. 

b. Transfer Function 

Let the input to the system be a plane wave, that is, g(t, r) = 
= r \ with frequency jo and spatial frequency vector u 0 . The output of 

the system is given by Eq. (1) as 



n no 

/2n/ t .<-,vr) e -j2*U,T-v 0 - ro) A(r>ro:r) dTdro! 



(3a) 



P(tr) = A(rro . r) dTdlo 

^ - ivT - /V\ 



(3b) 
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Fig. 2. Basic Geometry for the Interpretation of h{ r,r c :r). 



or 



VM = J 2 *V° t -' , ° ') S T .yS to .„{ h(T,r 0 ;r) } 



/=/o’ 

v=v 0 



(3c) 



where y represents the time-domain Fourier transform and _ v represents the 
spatial- domain Fourier transform. The subscripts under the symbol £ represent the 
variables involved in each transform. We define H(fv\ r), the system transfer 
function, as follows: 



H{fu\i) = h( T) r 0 ;r) } = J J < 



-j2w{fr-i/-T 0 ) 



^(r,r 0 ;r) dr dr 0 . (4) 



Using this definition, the output of the system can be written as 

ip(fi) = H(fo,v 0 ;r) ^ 2x (-4 ,_I/o ‘ r ) ; (5) 

the expected result from linear systems theory. 

The parameter r 0 in the integrations in Eqs. (1) and (4) represents the 
vectorial difference between the observation point r and the point source position r s . 
The integrations must be performed taking rasa constant, that is, by changing the 
source position r s so that r 0 assumes all possible values in R 3 , and the results, p(£,r) 
or H(fo,v Q \i), are valid for that "fixed" observation point r. 
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c. Response to a Time- Harmonic Point Source 

When evaluating the transfer function, it is convenient to use a 
time -harmonic point source as the input, <7(2, r) = 6(r-r s ) with frequency / c 

and location r s . The output, from Eqs. (1) and (4), is given by 



<p(t, r) = f f <5( r -r s -r 0 ) A(r,r 0 ;r) dr dr 0 , (6a) 

* - - /w 



<p(ti) = e 3 ^ T ^> r h( r,r-r s ;r) dr, 



(6b) 



<p(t, r) = ^ 27r ^' / 5 r _/{ h(r,r 0 ;r) } 



/ = /o ’ 

r c = r-r s 



(6c) 



or 



f{U r) = e? 2 } . (6d) 

Iq -IS 

d. Transform Relationships 

(1) F'ourier Transform Definitions and Notation. The time -domain 
Fourier transform of a function g(tj) = g(lx, y,z) and its inverse are given, 
respectively, by 



G(J,t) = d t _j{ g{t, r) } = J g{t, r) f 3 ^ l dt , 



(7) 
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and 



J oo 
-00 



( 8 ) 



The spatial Fourier transform and its inverse are given, respectively, by 



G(t,v) = { g(t, r) } = J g(t, i) (P^^dx , 

* 00 



( 9 ) 



and 



g(t.i) = G(t,i/) } 




e-j' 2 ™- r dv. 



(10) 



where dr = dx dy dz , dv = df* df y djz , and v-x = xjx + yfy + zjz . The above 
transforms involve triple integrals and can be split into three spatial transforms in x 
(i;)i V (/’) and z (jC). For example, the transforms in x (^) are written as 



G(tJ x ,y,z) = d T _^ { g(t,x,y,z) } 




(ID 



and 



g(t,x,y,z) = d' T [j y { G(ts,y,z ) } 



J oo 

e" i2x/ ^ . 



(12) 
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The transforms in y (fy) and z (j£) follow the same notation. The full space -time 
transforms can be written as 



G (U) = { 9 ^ T) > = d V-fy ^ > ' ( 13a ) 



or 



G(M 



n oo 

-m fVS 



2x(/t-i/-r) 



g(t,r) dt dr , 



(13b) 



and 



9{t-. r) = v { G{U) } = d}[ jd' T [j x d y [^ ^;ljr { G(f,v) } , (14a) 



or 



K , r! = JP^- r) o(MdJd». 



(14b) 



The above relationships are valid for the impulse response/transfer function pair if t 
is replaced by rand r is replaced by r 0 = (xb,J/o,2o), that is, 



Wf/r 0 ;r) = $ T -f{ h ( T »r 0 ;r) } , 


(15) 


H{t,v,i) = $ ro _ i/ { h ( r,r 0 ;r) } , 


(16) 
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W((, jft,!fo,Zo;r) = 5^ , { h(j-,xb,sb,%;r) } , 



(17) 



and 



= d T ,jd Io _ l/ { h(r,r 0 ;r) } . (4) 

The same applies for the respective inverse transforms. 

(2) Input -Output Relationships. During our developments, we will 

need to express the output of the system, £>(t,r), in terms of the transformed input 
function and the transfer function. From the Fourier transform definition 



9(t, r) 







(18) 



follows the transform property 



g(t-T,T-r) 




2x(f t.-v-i) 



c -j 2*^-^To) G(>) Jfiy 



(19) 



Substituting this expression into Eq. (1) yields 

00 

J'lxtft-v t ) e -j2ir(fr-i/-T 0 ) 

00 

x A(r,r 0 ;r) dfdv dr dr 0 , (20a) 
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p(J>r) = f l/ ' r ^G(/i/) f r#T,r 0 ;r) dr dr 0 df dv, ( 20b) 

^“OC^-OC ^“00^“ 

M 7, y; r) ' 



or 



pu r ) 







(20c) 



which is the desired relationship. 

Analogous to a time-variant system, for which the frequencies at 
the output are different from the input, a space-variant system causes the output 
spatial frequency vector to be different from the input. We will use /? to represent 
the output spatial frequency vector. The input -output relationship in the 
transformed form can be obtained by substituting Eq. (20c) into the Fourier 
transform definition 



Hf.P) = 




{ p(t,r) } e J l & T, dt dr . 



( 21 ) 



Doing so yields 



$(//?) = j j j v df] dv e 3 ^ & ^dt dr, (22a) 



■ OC “ CC “0C “ cc 



<!>(//?) = f f f G{rj,v) r) e ^ r j f e drdj/, (22b) 

" - <xr - vr - & W-oo / 



n 



= [ f | f G(tj,v) r) e ' j2 ^ l/ '^' r 5(/-f?) dr] J dr di/, (22c) 

" - OC* 7 - 00 W -OC ' 



4(M= J (22d) 

“ no 



2. The Acoustic Channel 

As long as we deal with small -amplitude acoustic signals, the governing 
wave equation is linear. We will restrict ourselves to time -invariant problems, 
where the properties of the medium do not change with time, and the source, 
receiver and the boundaries are motionless. 

Figure 3 is a pictorial representation of the acoustic channel. A source or 
transmitter rectangular array "couples" an electrical signal to the medium and a 
receiver array transforms the acoustic signal back to electrical form. Both arrays 
are assumed to act as linear filters. We will not be concerned here with the 
generation of the electrical transmitted signal, assumed given, nor with the 
processing done on the received electrical signals at each hydrophone position. The 
acoustic channel can be represented by the block diagram shown in Fig. 1, which 
can be broken down into the three basic components sketched in Fig. 3: transmitter 
array, medium and receiver array. Figure 4 shows these components in block 
diagram form. [Ref. 5;pp. 3-4] 
a. The medium 

The propagation of the acoustic signal through the medium is governed 
by the wave equation 
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Fig. 3. Geometry of the Acoustic Channel Problem. 
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Fig. 4. Block-Diagram Representation of the Acoustic Channel. 



1 a 2 

V 2 y?( «,r) - tp( t,r) = g( tj ) (23) 

<C(r) dt~ 

and appropriate boundary conditions, where p(£,r) is the velocity potential 
(m~s ), g( t.i) is the source term, which represents the rate of mass injection into 
the space per unit mass or rate at which fluid volume is added per unit volume 
( s 1) [Ref. 5;p. 1). These two terms are represented in Fig. 4 as the output and 
input to the acoustic medium block, respectively. The solution of Eq. (23) is 
obtained from Eqs. (1) or (20c) as 



= I I g(t-T,i-T 0 ) h M (r,r 0 ;r) dr dr^ , 

* “OC * 7 



(24) 



or 



•pM 




GUr) 



(25) 
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The transform of the output <!>(//?) is a function of the spatial 
frequency vector /?, which is, in general, different from v, due to the space-variant 
nature of the medium. The physical meaning, in terms of ray acoustics, is that a 
non -homogeneous medium (H is a function of r) causes change in the direction of 
propagation of the sound. 

The input -output relationship in the transformed form can be obtained 
from Eq. (22d), resulting in 







G(f.v) r) 

<x 



e - j2x(j/-/?)-r 



dr dv . 



(26) 



If the system is space -invariant, an unbounded homogeneous medium 
for example, then is independent of r and Eq. (26) reduces to the expected 
expression for a tirne-invariant. space-invariant linear system: 



$U0) = J* G(f,v) ®' T dr dv , (27a) 



= j G(f,v) H m (/v) f>(v-0) dv , 



(27b) 






(27c) 



b. The Transmitter 



The transmitter array, the first block in Fig. 4, is characterized by a 
complex aperture function A T (/r) (Amp 1 s relating the electrical input x(f,r) 
(Amp) to the source distribution g(«,r) (s' 1 ) through the expression [Ref. 5:p. 30] 

G(/r) = X(f,r) A t Ut) • (28) 

3 -1 -1 * 

The far -field directivity function (m Amp s ) is related to the complex 

aperture function through the transform [Ref. 5:p. 38] 

VM = S r .,d -VM}- (29) 

Using Eq. (29), the spatial Fourier transform of Eq. (28) is given by 

G{f,v) = * D t (/ } v) , (30a) 



or 



(?(/*')= X(J,a) D T (f,i/-a) dct , (30b) 



where the symbol (*) denotes a three-dimensional convolution with respect to v. 
Equations (30a) and (30b) give the relationship between the electrical input and the 
source of the acoustic signal. The relationship between the input electrical signal 
and the velocity potential (the output of the medium block of Fig. 4) is obtained by 



1C 



substituting Eq. (30b) into Eqs. (25) or (26). Let us introduce the simplifying 
assumption that the electrical input signal is independent of position, that is, 
jr(f,r) = x(() and X(f,v) = X(f) 6(v) . This is justified by the fact that the 
complex weights applied to the electrical signal (array shading and beam steering) 
are taken into account in D^(f, v), the far -field directivity function. With that 
assumption, Eq. (30a) becomes 

G{U) = XU)6(»)lD T U») t (31a) 



or 



G(U) = XU) £> t ( j>), 



(31b) 



which, when substituted into Eq.'(25), gives the required relationship between the 
input electrical signal and the velocity potential at the output of the medium block, 



?{Lv) = j j X(f)D T (j,v)H M Uv\T) df. 



(32) 



A further restriction we impose to our problem is that the transmitter array is 
composed of point sources, has rectangular shape, is parallel to the XY plane, as 
shown in Fig. 3, and is centered at position r s = (:%, Jts.-Zs)- The array has 
dimensions Mt * Nt (first dimension along the X direction and second dimension 
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along Y), with interelement spacing d x t meters along X and dy t meters along the Y 
direction. Then, the far -field directivity function for Mt and Nt odd is given 
by [Ref. 5;ch. 4] 



M t -l N t -1 

— 5“ T 

D r (M = £ X W- f >* )2 " , ’ rpq - (33 > 

Aft-1 _ Nfl 



where r^ ^ = (pd x t + £s. qdyt + y s , •%) represents the position of each element and 
c (/) = <z pq (/) is the complex weight associated with the element at 

position r pq . 

c. The Receiver 

The receiver array, the last block in Fig. 4, is characterized bv a 

r 

complex aperture function A (f,i) (V sm °) relating the velocity potential <p( f,r) 

R 

-3 

to the electrical signal y( t, r) (V m °) through the expression [Ref. 5:p. 31] 

= HJ,t) \U, r). (34) 

As stated earlier, we are not concerned with the signal processing done 

on the received signal. Our purpose is to predict the waveform at the output of 

each receiver array element. So, we consider the term A (f. r) as pertaining to a 

single transducer element. The quantity y(i, r) is the electrical signal distribution 

-3 

(V m ) due to the velocity potential p(£, r) at a point r on the transducer. In order 
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to obtain the total electrical signal y(t) at the transducer terminals, we must 
integrate y(t,r) over the whole transducer volume, that is, integrate with respect to 

r, 



n/) = J"V«o *, 

“ oc 

Y(f)= f *{tT)AMv)dr. 



(35) 



ky 



For a point sensor located at r r , the complex aperture function is given 



A R (f,i) = A{f) 8(i-i r ) , (36) 

which results in an electrical signal at the terminals given [ see Eq. (35) ] by 

Y(f) = Hfcr) A(f). (37) 

If the transducer is a pressure sensor, then .4(/) is proportional to the frequency /, 
that is, and Y(f) ~ ~ pressure . For simplicity, we take A(f) 

as unity, and a pressure sensor can be simulated by inserting a suitable filter into 
the block diagram of Fig. 4. So, the output electrical signal will be numerically 
equal to the velocity potential. The receiver array, like the transmitter, is a 
rectangular array of point sensors of dimension M r * N r , parallel to the XY plane, 
with interelement spacings d xr and d yT meters, centered at r r = (ly, J/r> %)• The 
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output electrical signal at each element (m,n), given by Eq. (37) with A(f) = 1 
and T r = i mn = (m, dxr + Xr , n dy r + y r , Zr), is 

Vinn(/) = <J>(/r mn ), (38) 



or 



^mn ( ^ ^ r mn ) * 



(39) 



d. Overall Transfer Function 

The overall transfer function relates the output electrical signal J/ ran (t) 
at each point sensor ( m,n ) of the receiver array to the electrical input signal .r(t) at 
the transmitter array. Substituting Eq. (39) into Eq. (32) we obtain, for the output 
signal, 



«,„(<) = f f A'(/) D t (M Jf, (40a) 



J oo /• oc 

A(/)J D^U,v) (40b) 



■CC ”00 



The last expression can be rewritten as 



” 00 



(40c) 



20 



or 



sUO = K f {xu) «(A™)> 



(40d) 



where H(j,T mn ) is the overall transfer function, given by 



H^f,- !r m „) r-' 2 «' r "« ( fo 



(41) 



Note that Eq. (40c) is valid only when the input electrical signal is independent of 
position and the receiver is a point sensor with unit frequency response, that is, 
A(f) = 1 . On the other hand, Eq. (32) assumes only the independence of the input 
electrical signal with respect to position. 

B. THE RANGE INDEPENDENT ACOUSTIC CHANNEL 

In this section we will derive the expression for the overall transfer function for 
the particular case of a range independent medium. With reference to Fig. 3, the 
range independent medium characteristics change only with the Y (depth) 
direction, and the boundary conditions are independent of r and z , that is, the 
problem presents cylindrical symmetry. 

1. Medium Transfer Function 

In order to obtain the medium transfer function, we will solve the wave 
equation using a time -harmonic point source, that is, 

V>(t,r) - -7 <5(r-r s ) , (42) 

c~(y) dr 
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with suitable boundary conditions . The solution, as shown in Eq. (6d), is of the 
form p(J,r) = £(r) J w ith £(r) satisfying the Helmholtz wave equation 

V 2 f(r) + t 2 (j,){(r)=«(r-r«), (43) 



where k(y) = 2irf /c(y) is the depth dependent wave number. Applying the spatial 
Fourier transform with respect to rand z to Eq. (43), we obtain 

k-Au) S«,s a) + 4jE(jUjO = + w 6(y-y,) , (44) 

dy 

.0 . . .o . . o , . o .0. _ . 

where ky(y) = k~(y)-4 ;r"(jk“-i- j£~) and is the spatial Fourier transform 

of £(r) with respect to i and 2 . Solutions to Eq. (44) are usually obtained by 
approximate methods, most notably the WKB approximation. We can write the 
solution to Eq. (44), without loss of generality, as 



^(ic ,y,i) = a(y) < A t 



■A ( y) + B JOyiy)) J2 j(J k Xs + f z Zs) 






(45) 



-S/'* 

4 > 



The (transformed) velocity potential, obtained from Eq. (45) and the assumed 
fUr) = f(r) is given by 

+( = H(A,!/,0 . (46) 



OO 



The velocity potential <p(i,r) is given by the inverse Fourier transform of Eq. (46) 
with respect to rand z , namely 






Lx -j'ljtz . r j C 

Jx e J J * d& dh . 



(47) 



Next, we will relate these results to those obtained in Section II-A-lc, the 
response of a linear time-invariant, space-variant system to a time-harmonic point 
source, as given by Eq. (6d), which can be written as 










(48) 



where r 0 = (jb, J/o, A.) = (x-j%, y-y S) z-z % ) and dv = dj x df y df 2 . In order to 

compare Eqs. (47) and (48), we rewrite Eq. (47) using the term 'h defined in 
Eq. (45) as follows: 



7 (t.T) = j2 * lJ fr* ■ JMIxXs + /z2s) e~ j2rI * Z df x df z . (49a) 



<p(t, r) = 





2x/xTc e' j2 '^dj x dj z 



(49b) 



Comparing the Fourier transform expressions in Eqs. (48) and (49b). the term 
obtained from the solution of the inhomogeneous wave equation when the source is a 
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unit amplitude, time-harmonic point source, is seen to be related to the medium 
transfer function through the transform 



y(f,kyo,&y) = r y ' o _ fy { }, 



(50a) 



or 




" 00 



(50b) 



where y has been substituted for r as an argument of H ^ in order to stress the 



depth -only dependence of the medium transfer function. Notice also that the 
arguments of have been written in accordance with the linear systems notation 



introduced in Section II -A - 1 . 

The transfer function can then be obtained as the direct spatial Fourier 
transform of 'fM J-fx-yo.fcy) with respect to y o — y-y s - In the next section we will 
see that this last step is not necessary, so long as the transmitter array is composed 
of point sources. 

2. Output Electrical Signal 

The overall transfer function for a range independent medium can be 
computed from Eq. (41) by substitution of the expressions for D t (/i/) and the 
medium transfer function. As we just mentioned, it will not be necessary to 
compute the medium transfer function if the transmitter elements axe point sources. 
The directivity function for the transmitter array is a sum of terms of the form 
Cpq(f) & ~ pq j where r pq = (z p , y q , is the position of each transmitter 
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element. Substituting that general summation term into Eq. (41) and expanding 
the exponential factors yields the result that the overall transfer function is a 
summation — over p and q — of terms of the form 



'pq 



U )j~d Hvtvw H t u,u,ky°) , 



or 



C pq(/) 



I 

1 



rjhfJx. 



T 

m V e 



■i 2<r/ y (!/ n -!/ q ) d „ . 



where x m = x r +m d xl and y n = y r +n d yr are the (x,y) coordinates of the receiver 
array elements, and the (xy) coordinates of the transmitter array elements are 
= W d xt and y q = y^+q d yt . Recall from Section II -A -1 that the parameter 
r 0 is the vectorial difference between observation point and source position. The 
integral in this last expression is the inverse spatial Fourier transform of the transfer 
function //(/£, fv,fi'-yn) with respect to jx , fy and j£ which yields a function of 
(xo: Vc- -^c). where Xq = , y 0 = y n -y q ■ and Zo = • But the transform 

with respect to Jy is already given by Eq. (50b) as the function $ (JJx,yo,&y) 
obtained from the solution of the wave equation. As a result, we can rewrite the 
above last expression as 



C pq(/) 



n OO * • 

^UJx.l/o -ky) e e 3 ~'*h~odfy L dj, z . 
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Therefore, the overall transfer function for the range-independent medium is given 

bv 



ft-1 ^t-1 

T" ~~T 



i i ^(/)fr 

#t-i -<» 

P = --y = 2~ 



VU&yot&yn) e~ j2 *k X o 



X e -j-*& Z o df x dk . (51) 



The output electrical signal at each receiver position is given by Eq. (40c), 
repeated here for convenience: 



y„(l) = f «/) W(/r, M ) 



(40c) 



where H( f, r^) is given by Eq. (51). Again, we stress that Eq. (40c) incorporates 
three assumptions regarding the linear system: 

• The input electrical signal is independent of position r, that is, x(i,r) = x(t) . 

• The receiver array elements are point sensors, that, is, A (f, r) = A(f) £(r-r r ) . 

• The sensors have unit frequency response, that is, A(f) = 1 . 

Equation (51) assumes, in addition: 

• The transmitter array elements are point sources. 

• The medium has range independent characteristics and boundary conditions. 

• Transmitter array geometry is as described in Section I -A -2b. 



a. An Alternative Representation 

We will now compare our results, as given by Eqs. (51) and (40c), with 
the results obtained by solving the wave equation for a range -independent medium 
in cylindrical coordinates and written in terms of a Fourier-Bessel integral. 
DiNapoli and Deavenport [Ref. 4:p. 80], for example, use that representation as a 
starting point for the description of numerical methods applied to range independent 
problems. 

The integral in Eq. (51) can be seen as the overall transfer function due 
to a single point source. From this point of view, H(f, r ran ) is just the weighted sum 
of the transfer functions due to each transmitter array element, as one would expect 
for a parallel connection of linear systems, These elementary overall transfer 
functions H pq ( fj mn ) are given by 





OC 

OC 



:-j'^k X 0 r -j~xfz Z O 



dk dk 



(52) 



The overall transfer function can then be written as 



A/t-1 
— T~ 

W r nJ = X 

M t -1 
V = 2 “ 




C pq(/) W pq(/ r rcJ • 

Aft- 1 
~T 



(53) 



Let us make a change of variables in the integral in Eq. (52), exploiting 
the fact that the dependence of ^{Ik^kVn) with k and k occurs due to the 
factor ky(y) = £ 2 (j/)-4x 2 (/ x ^+ $) in Eq. (44), that is, V is a function of (jk~-r k*)- 



Introduce the variables / r and 7 such that 



/x = X coy 7 and f z = f T sin 7 , 



(54) 



which is a change from "rectangular" coordinates (X,X) to the "polar" coordinates 
(X, 7)- With this change of variables, Eq. (52) becomes 



flp, = r*(«,w*> f, cos 1 + si " **, a , , 

J ^ fi 



(55) 



where 



^{U v y o' ; yn) = y(fJx,yo,kyu) 



X=X cos 7 
X=X 5l ‘ ra 7 



(56) 



Equation (55) can be written as 



f/ pq(/ r :rJ = 



2 j 

[ / r I 4 , 

J (\ J f) 



(57) 



where sin a = x 0 / R 0 , cos a = z 0 / R 0 , and R 0 = \J x^ + . Because of the 

periodicity of the sine function, the term a can be dropped from the exponent 
without altering the result. The last integral can then be recognized as the 
zero-order Bessel function of the first kind [Ref. 6:p. 184], that is, 
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(58) 



JM 



lf 2 '. 
= ^\ e 



-jz sin 



7 <*7, 



which leads us to the Fourier -Bessel integral representation 



= 2x f 



= 2' nu,y°;y ») 4(2»/A) /, <0f ■ 



(59) 



Equation (59) is the same as Eq. (3.1) in [Ref. 4] if we perform the 
trivial change of variables k r = 2x£. There, T is interpreted as the depth dependent 
Green's function. Equation (59) is the starting point to several approaches to the 
range independent problem, both numerical and analytical [Ref. 4:pp. 80-134]. 
Using the Fourier -Bessel integral representation, the overall transfer function 
H(f r rin ) can be written, from Eqs. (53) and (59). as 

M t -1 A’ t -1 

T 

-t. lu - ^ 

H(/r„„)= £ £ c pq(/) ~ J j ( 60 ) 

M t -1 iVfl 

* = tf = "-T- 



C. THE LAYERED WAVEGUIDE 

In this section we will derive the function T for the case of a three fluid medium 
waveguide, commonly referred to as the Pekeris waveguide. 
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1. Statement of The Problem 

Figure 5 shows the relevant physical aspects for our derivation. The point 
source is at position ( x s ,y s ,z s ) in medium II and the receiver is at (x,y,z) in the same 
medium. Medium II, with characteristic impedance p^c^, is bounded at y = 0 , the 
surface, by a fluid with characteristic impedance />,c, and at y = D , the bottom, by 
another fluid with characteristic impedance /> 3 c 3 . As described in Section B-l, in 
order to obtain 4q we need to solve the Helmholtz wave equation, Eq. (44), and 
isolate the terms with y dependence from the solution, as shown in Eq. (45). An 
inspection of those two expressions show us that by setting x s and z s to zero, the 
solution of the wave equation will be automatically the function The wave 
equations that describe the propagation in the three media are: 



iyi 4*, 



dy“ 






= 0, 



( 61 ) 



o 

Ay 2 4* 



i 2 

o H 77 

if 



^2 = t(y-ys) , 



(62) 



and 



o 




^3 + 7^2 ^3 = 0 > 

dy 



(63) 



2 2 2 2 2 

where k yi = }^) , k i = 2 and the subscripts 1, 2, and 3 refer to 

the three different media. The argument of 4> has been dropped for simplicity. We 
seek the solution 4' 2 because we shall be concerned only with those problems where 
both source and receiver are in medium II. The boundary conditions are the 
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NJ 





Fig. 5. The Three Layer Waveguide. 
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continuity of the acoustic pressure and the normal component of the acoustic 
particle velocity at each boundary. Acoustic pressure is related to the velocity 
potential by 



Pi(/) = i2*/ft*i, (64) 

and the normal component of the acoustic particle velocity at the boundaries is 
related to velocity potential by 



U = — \]> ■ 

^ m , 1 



dy 



( 65 ) 



2. Solution 

With the source located in medium II, the solutions in media I and III must 
be in the form of propagating waves moving outward from the boundaries. In 
media II, there are waves propagating in both directions — increasing and decreasing 
y — due to the reflections at the boundaries. The expected form of the solutions are 



2a - ^2a 



o 

\/l 

H 

& 


(66! 


<■ A'= J - B 2 , y , 0 < y < y t , 


(67) 


A’ » + B 2b e'*y^ y ,y,<y<D, 


(68) 
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and 



* 3 = B 3 e-fy* y ,y>D. 



(69) 



The solution in medium II has been split in two in order to account for the point 



source, as usual in the Green's Function determination. Two more boundary 
conditions are necessary: 

• Continuity of the solution \k 2 at y — y s . 

• Discontinuity of the first derivative of 'I' 2 at y = y s . 

The value of the discontinuity can be found by integrating Eq. (62) over a region 
I U~y s I < f and taking the limit as e -* 0 . Performing this operation and taking into 
account that, in the limit, the integration of the continuous function is zero, we 
obtain 




(70a) 



or 




(70b) 
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a. Boundary Condition at y = 0 

In order to obtain the expression that represents the continuity of 
pressure at y = 0, we equate the value of pressures due to and 'f' 2 a - Substituting 
the assumed solutions — Eqs. (66) and (67) — into Eq. (64) yields 

P\ = P‘1 (^2a + ^2a) • (71) 

For the normal component of particle velocity, Eq. (66) and the assumed solutions 
yield 



M 



.u 



= X v2 



'.^2a"^2a^ 



(72) 



Substituting Eq. (71) into (72). is eliminated and the relationship between the 
amplitudes of the solution 'K a is found to be 




P i 

— k - k 

p 2 V S i 

- k + k 
p 2 y- + *y» 



(73) 



When L is real, corresponding to propagating waves in medium II, the factor /? 21 is 
physically interpreted as the plane wave reflection coefficient at the surface. 

b. Boundary Condition at y = D 

Following the steps of the previous section, using the forms of the 
assumed solutions for ^ and 4s b at y = D , we find 
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Pi A, e-*>* D = Pi (A* + Bj b ) , 



(74) 



and 



-* y3 = l y2 (/) 2b e'^ D ) . (75) 

Substituting Eq. (74) into (75), .4 3 is eliminated and the relationship between the 
other two amplitude factors is found to be 



“2b _ ,• 

= FUo e y 

% 



2k. .^D 

y z = 



* * . k 

P, V )'3 

/J . 

~ k + k 

Pv K y- j 3 



e -j2k y2 D 



(76) 



When is real, has the physical meaning of the plane wave reflection 
coefficient at the bottom. 

c. Source Condition 

Using the results of the two last sections. Eqs. (73) and (76). the 
solutions in medium II become 

ty 2a = /U a (A'*y + /i,, e ), 0 < y < y s , (77) 



and 



^2b — ^2 



(tf 23 e'j- k y'- D ei k yi y + e"^y2 y ), y s < y< D 



(7S) 
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The condition of continuity of 4^ at y = y s , applied to the above two equations, 
gives the relation between the amplitudes A 2a and B^: 



A 2k r> D J k y* y ° + e-i k y*y° 

v, + Rn C -J* y2 Vs 



(79) 



Finally, apply the condition on the derivative of 'I' 2 at y = y s , as given by 
Eq. (70b), to obtain: 



^ „ - j 9 jfc „ D ik o 

^2b 1*23 C • y " T2_ ^ Y2 



is - f ”$y2 is 



M 2a (A^./i 21 ) = i_ 






( 80 ) 



V2 



which, together with Eq. (79), forms the necessary system of equations to obtain the 
two amplitude factors. B 2 1> and .4 2a . Upon substitution of the expressions for B 2b 
and ,4 2a into Eqs. (77) and (78), we obtain the solutions ty 2a and T 2b , which form 
the transformed transfer function T (J,k, ycihw) we seek. Written in an 
abbreviated form, the final expression is 

R a e- j2k y°- D e jk v'- y > + c-^ y> 

^UJx,yo-ky) = j— 1 

- k r- 1 - Ft,, R r , fi - k r'- D 

( 81 ) 

where 0 < y < D , 0 < y s < D , and, as usual, y> is the greater and y< the smaller 
between y and y s . 



3. Special Cases 



In order to check the solution for the layered waveguide, we will now 
derive the expression for the output electrical signal in two special cases whose 
analytical solutions are known. These results will be used to check the computer 
implementation. In both cases, the input will be a time-harmonic point source. 
a. Unbounded Homogeneous Medium ( Free Space ) 

An unbounded homogeneous medium can be represented by Fig. 5 if 
the sound speeds Cj and densities pj are equal for the three media. In this case, the 
reflection coefficients, as given by Eqs. (73) and (76), become zero. The 
transformed transfer function — Eq. (81) — reduces to 

*(LLyo<kv) = 3— . (82a) 

)b 

*w/ty2 

^(fkyo-kv) = j— e' jk y^y> ' y <] . ( 82 b) 

■jk y2 

or 

*(f,k!h,ky) = j— t~ jky ^ y ' ^ sl = j— e' jk y^ . ( 82 c) 

2Ay, ~ k y2 

Note that the right hand side of Eq. (82c) is expressed in terms of y Q only, meaning 
that the unbounded homogeneous medium is space -invariant, as expected. 

The overall transfer function H(Jj) for the case of a point source is 
given by Eq. (51) with M t = N t = c,,(/)= 1. Substituting Eq. (82c) into Eq. (51) 
yields 



HUt) = f ( “j -i- it . ( 83 ) 

" -~y- oo 2 ^y2 



The above integral is the expansion of a spherical wave into plane waves [Ref. 7:p. 
B-5] and the resultant overall transfer function is 



HU, r) = - 



- jk 2 R 
4 jR 



(84) 



f~n T 7 / 75 75 7 

where R = || r-r s || = V -r £ -r = \/( 

The output electrical signal ?/(/) — or velocity potential <p(t, r), in our 
case — is given by Eq. (40c). Substituting Eq. (84) and X(/) = 6(f-jo ) — for the 
time-harmonic source — into Eq. (40c) yields 



y(t) 



_ . _L_ e i’ 27 r X-, / - KR ) ; 
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(85) 



where Ay = /c 2 . liquation (85) is just the expression for a spherical wave, since 

R is the distance between source and receiver. 
b. Surface Reflection 

Consider the case where media II and III are the same, that is, c 2 = c 3 
and p 2 = Pi ■ As a result, is zero and reduces to 



^U,kyo,ky) 



2k 



y 2 



e -i^y> + /?2| e -i^y< > 



(86a) 



•aUkvo.ky) = /— (t' lk r> [y> ' y<) + flj, + »<>), (86b) 

2t y2 

or 

1 '(},kyo,ky) = y-i— ^ + flj, + ) , 0 < y, 0<J/ S . (86c) 

2* yJ 

Comparing Eqs. (82c) and (86c). we see that this problem can be treated by 
considering a homogeneous medium with two point sources, one located at y . with 
unit strength, and the second located at -y s with strength FU V The output 
electrical signal due to the first term in Eq. (86c) has already been computed in the 
previous section, and is given by Eq. (85). In the second term, we have the factor 
FLy |. also a function of f K and J z through the dependence on k yl and k y2 . When 
performing the inverse spatial Fourier transform to obtain H(j, r), it would be easier 
if R 2l were constant. If we consider the boundary I -II as pressure release , that is. 
pjp 2 = 0, then, by Eq. (73), /v,, = -1. In this case, Eq. (86c) reduces to 

y(f,jx,yo,.L;y) = j — — (e’A^W-e ■*W y ol) . (87) 

O 

- * v 2 

where y ' = y-{ -y s ) . From the results in the previous section, the output electrical 
signal is, by superposition. 



?(0 = 



_L 

UR UR' 



( 88 ) 
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7 s — 7 
+ Vo + Z o 



! ^ — 5 2 

= v( x-x s y+(y+y s ) +(z-z s ) . If the distance 



where R' = || r-r^ || = 

R is large compared to the depth of the source y s , a simplification is possible. From 
Fig. 6, for large R , Rz R S -AR , R' a R s + A R and A R z y s sinQ = y s y r / R s , 
where R s is the distance from the receiver to the midpoint between the source and 
its image. With these approximations and letting R z R' z R s in the amplitude 
factor, Eq. (88) becomes [Ref. 7:pp. 170, 408] 



y{ 0 



= - w 



sm 



2 x/L 



( k o y_s Vy \ 

R. 



(89) 



D. SUMMARY 

The transformed transfer function yo,jz',y) for a range-independent 

medium is the solution to the inhomogeneous Helmholtz wave equation 



o 

ky(y) V + # = 6(y-y s ) 

dy~ 



(90) 



2 2 2 2 2 

with suitable boundary conditions, where ky = /7-4x (jk + h") and &=2x//c. 

Under the conditions stated in Section II -B -2, related to Eq. (51), the overall 
transfer function for a range-independent medium is given by 



M t-1 
T 



N t -l 
T 



H U> T m n) = ^ f * U>J*>1to>kyn) 

M t - 1 _ N t -1 '» ■» 

” " ~~T 

x t -3-*S z z o df x djz , (51) 



V" 9 = 
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Fig. 6. Geometry of Surface Reflection. 
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where the transmitter array is rectangular, centered at (x s , j/ s ,z s ), parallel to the XY 
plane, with odd dimensions Mt x Nt , complex weights c pq (f), and with 
elements positioned at x p = x s +px , y q = y s +qy , and z = z s . Also, x Q = x m -x p , 
y 0 = y n -y q , and z 0 = z T -z s . The overall transfer function H(f, r m ) can also be 
expressed as 



The output electrical signal y^it), under the same conditions stated for the 
overall transfer function, is given by 



A/t-1 



Nt-1 




* J 0 (2wf r R 0 ) f r df r , (60) 



where 




^ (ftlrtVoWnl — 'K (f>jx,yo,jz] 2/n) 



/x=/ r 

fz=fr sin l 



(56) 




oc 



(40c) 



where H(f, r mn ) is given by Eqs. (51) or (60). 
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The layered waveguide illustrated in Fig. 5 is characterized by the transformed 
medium transfer function 



a,, y> + e -i k s>y > 

y{lk,yo,ky) = j— : 

2 *ys l- , ^23 2ky2 ° 

+ (81) 



which is valid for 0 < v < D , 0 < < D . and where R n , is eiven bv Ea. 173) and FL-, 

by Eq. (76). If all three media have the same acoustic characteristics, then Eq. (81) 
reduces to 



'Hik-yU-v) = i— , 

O 1 

4s fV.,0 



($'2c) 



and the corresponding output electrical signal is given by 

*(()=-- (65) 
4 wR 

where R is the distance from source to receiver. In the case where media II and III 
are the same and the surface is pressure release, the function ^ is given by 

= i— (e^yM.e-j^y'ol) t (87) 

2k y2 
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where y' Q = y + y s . The output electrical signal, when the source -receiver distance 
is large compared to the source depth, is given by 



y(t) = 



-,-_L (l!ii 

2 if R, { R, 



(89) 



where R s is the distance from the receiver to the midpoint between the source and 
its surface image. 
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III. IMPLEMENTATION AND RESULTS 



In this chapter, a computer implementation of the range independent medium 
equations is described. Although only the medium transfer function of the layered 
waveguide is implemented, the modular nature of the linear systems approach 
enables us to analyze the results and draw conclusions valid for the more general 
case. Equation (60) illustrates this point. The summations and the coefficients 
c pq(/) * n ^hat e q ua -tion represent the transmitter array. The Bessel function 
J 0 ('2~j t R 0 ) is the range dependent factor. The function ^ (f,f^,y 0 J z ^y ) — from now 
on referred to as the medium transfer function, for simplicity — is the (range 
independent ) medium dependent factor. 

Our purpose is to have a working algorithm to test the results from the previous 
chapter. Care has been taken to enable the implementation of other medium 
transfer functions, allowing for depth dependent speed of sound, so that the program 
can be used as an academic tool. Speed and style were secondary objectives. 

The implementation is based on Eq. (60), the Fourier -Bessel integral 
representation. 

A. IMPLEMENTED EQUATIONS 

1. Medium Transfer Function 

The medium transfer function 'I' (// r ,y 0 ; y) is an oscillatory or exponential 

2 2 

function of k y2 — as seen in Eq. (81) — which, in turn, is a function of (f x + f z ) = 
= j~. The function 'I' is easier analyzed in terms of k y2 than in terms of / r = 
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= J ifi^y r - (k^2x? , which is nonlinear. Therefore, if will be treated as a 
function of k y2 , and the integral term in the expression for the overall transfer 
function, Eq. (60), will be modified accordingly. 

a. if as a Function ojk^ 2 

The function is already expressed in terms of k y2 in Eq. (81). The 
reflection coefficients, on the other hand, are expressed in terms of k yl , k y2 , and k yS . 
In this section we establish the relationships between k y2 and the quantities / r , k yl , 
and k y 3 . These relationships will enable us to represent the reflection coefficients in 
terms of the variable k y2 , as well as to change the variable of integration in the 
expression for the overall transfer function to k y2 . 

The correspondence between k y2 and f r is one-to-one. When j f < (//c. 2 ) , 



r — 5 — 5 

k y2 is taken as the positive root c 2j\J (/ fc^Y- f r , a consequence of our choice of 
solutions, Eqs. (67) and (68). When / r > (//eg), k y2 is imaginary, corresponding to 
an exponentially decaying ty. When all factors in Eq. (81) are multiplied out, is 
seen to be composed of a sum of factors of the form 



A e~ J ' k y* a 



) 



where a is a nonnegative real number. In order to have ^ be a decreasing 
exponential, the negative root must be chosen when k y2 is imaginary. Similar 
reasoning yields the same results regarding £ yl and A y3 , which can be summarized as 
follows: 
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*yi = 



^ ~ ^ , *r < 



(91) 



-J 



*r ~ *i , * r > 



where i — 1, 2 or 3 , A r = 2 and ^ = 2wf fc^. Using Eq. (91), the wavenumber 
y-components £ y , and i y3 can be expressed in terms of k y2 as 



*yp ~ 



+ 






~ + k y 2 , k T < kp 



a! P 



2 o o 

"■2 — 2 ~ l - > 



(92) 



I — 2 2~ 

where p = 1 or 3. and k { = ^ k 2 — k y2 . Using Eq. (91), the integrations over / r 
can be transformed into integrations over k yi . Using Eq. (92), the reflection 
coefficients /? 2( and R n , as given by Eqs. (73) and (76), can be expressed in terms 

of Ayo. 

b. Air, Water, Fast Bottom Waveguide 

The present implementation of the medium transfer function assumes 
that Cj< Co. When the surface is almost pressure release, which is the case of the 
air-water interface, and the bottom is fast, c 3 > Gj , the medium transfer function is 
composed of many sharp peaks, which presents some difficulties for integration 
routines. 

(1) Poles of the Medium Transfer Function. When the surface is 
pressure release, that is, R 2i = -1 , Eq. (81) reduces to 
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VU,k y2 ,yo]y) = 



ft,, c -j- k y2DJk y ,S> + e -jk s2 y> 

sm( k 2 y< ) , 

; V 1+ fl !s e' j2i r2 D 



(93) 



which is valid for 0 < y < £> and 0 < y s < D . If the bottom is fast, k^ > k 3 , and, 
from Eq. (91), there is a range of values of k r for which k y2 is real while & y3 is 



imaginary, that is, k% < k T < or 0 < 
of k y2 , the reflection coefficient 3 , given by 




For this range of values 



^23 — 



p 2 k _ k 

p<2 S 2 S' 3 

~ L r 

^ *y 2 + Ss 



(94) 



becomes a unit magnitude complex number, 



R« = 



(95) 



where 



ft 7m { *ys ) 

td 71/ (jj — ~ ' ——————— 

ft S 2 



Pi 

ft 




(96) 



In this case, the denominator of the expression for ty, Eq. (93), may have zeros, 
given by 



1 + ft,, = 0 , 



(97a) 
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1 + e ~j 2 (k y2 D ' 0 ) . 0 . 



(97b) 



Therefore, 



<p - k 2 D = (2n + 1) - , n = 0, 1, 2,... 

y 2 



(98a) 



or 



tan k y2 D 



-1 



^3 



y2 



tan 4> p 2 



ft . ft . ft 

*2 r 'y2 h 'S 



(98b) 



lirVi irli rn r> kn nvnrnor 

vv ii iv, j i ^<xii uc CA|'iroocu oo 



Pz 0 

tan 9 = 



P 2 



, « 2 -<P 



(99) 



■ ; v 

where 9 = £ y2 D and a = D ^ k 2 - fcg • Solutions to the above transcendental 
equation in terms of k y2 = 9/ D are the poles of the medium transfer function, for 
the case of a pressure release surface, Eq. (93), and a fast bottom. The physical 
interpretation of these poles is that they represent the discrete eigenvalues k y2n 
corresponding to the trapped or normal modes in the waveguide[Ref. 8:pp. 430-440]. 
The graphical solution of Eq. (99) [Ref. 8:p. 438] reveals that: 
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• Poles exist only if a > t/ 2. If a > x , the poles will be separated by A & y2 » ar / D. 

• The total number of poles is given by the integer number equal to [(a/r)- 0.5] 
or, if the result is not integer, its nearest larger integer. 

When the surface is almost pressure release, that is, R 2l « -1 , 

the medium transfer function presents peaks at values of k y2 corresponding to the 
poles computed as above. 

2. Overall Transfer Function 

In order to change the variable of integration in Eq. (60), the interval of 
integration f r € ( 0 , oo ) is subdivided into ( 0 , J / c^) and ( / /cj, oo ). After the 
change of variables according to Eq. (91) and a rearrangement, Eq. (60) becomes 



u( f, 



mn 



_r i 

Jr 



V V 
P 9 



- i n 
pq KJ ) 



\T< ( ( U v ■ ,. ^ T t L D \ L 
v VJ) A y 2'Po’^n^ J o\ ,v r il o> K y2 



dk 



V2 



n nn\ 



where the path of integration T consists of the imaginary axis over the interval 
(-jo c , j 0) plus the interval (0 , f / c) on the real axis. The summations over p and 
q are as indicated in Eq. (60). 



B. THE SUBROUTINES 

In this section we describe the subroutines and present the results for some 
examples. Also, a description of the main program and some auxiliary subroutines 
is given, which will help to clarify some aspects of the implementation 1 . 



•The main program and auxiliaxy routines were written by Dr. Lawrence J. 
Ziomek, at the Naval Postgraduate School, as part of an ongoing research project in 
model -based signal processing. 
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1. Main Program and Auxiliary Routines 

The main program's simplified pseudo-code is shown in Fig. 7. It consists 
basically of calls to subroutines. The readin subroutine reads the input data, 
consisting of transmitted signal characteristics, transmitter and receiver array data, 
medium characteristics, and control data, including the ocean medium transfer 
function (OCNTF) variable used to select a transfer function, when more than one 
type of transfer function is implemented. The ccq subroutine generates the input 
electrical signal x(i). The chfmn subroutine computes the overall transfer function, 
Eq. (100), and its estimated error. The phfmn subroutine provides a tabular output 
for the values of the overall transfer function and its error. The calyce subroutine 
computes the output electrical signal, according to Eq. (40c). 



Main Program /computes output waveform of an acoustic channel/ 
EXECUTE Readin /read input data/ 

EXECUTE Ccq /generate input electrical signal/ 

EXECUTE Chfmn /compute overall transfer function/ 

IF (flag print is TRUE) THEN 

EXECUTE Phfmn /print overall transfer function/ 

END IF 

EXECUTE Calyce /compute output electrical signal/ 

End Main Program 



Fig. 7. Simplified Pseudo-Code of the Main Program. 
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a. Subroutine Ccq 

The subroutine ccq generates the input electrical signal i(t), using the 
complex envelope representation. From the specifications read as input data, this 
subroutine computes the time samples of a reference complex envelope and its 
Fourier series coefficients. The Fourier series is truncated, yielding the complex 
envelope of the transmitted signal. 

(1) Complex Signal Representation. An amplitude and angle 
modulated carrier r(t ) can be represented by its complex envelope f(i) as [Ref. 

5: PP . 177, 182] 



r(t) = Re{ f(t) J } , 



( 101 ) 



where 



V) = «(0 <P l) , 



(102) 



a(t ) is the (real) amplitude modulating signal (Amp), 9(t) is the angle modulating 
signal (rad), and f r is the carrier frequency (Hz). The Fourier transform of the 
signal is related to its complex envelope transform by [Ref 5:p. 187] 



R(l) = 04 RU- 1 ) * «*(-/- /c)l • 



(103) 



For a rectangular -envelope LFM pulse, a(<) and #(<) are given by 



a(l) = 



'A , |t | < Tp/2 
.0 ,|1|> Tp/2 



(104) 
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and 



3T 5W 0 
6(t) = t 

T? 



Dr ?, 



(105) 



where A is the (constant) amplitude (Amp), Tp is the pulse length (s), 5V is the 

-9 

swept bandwidth (Hz) and Dp is the phase -deviation constant (rad s ). Both 5V 
and Dp are allowed to be either positive or negative, corresponding to the up chirp 
and down chirp LFM pulses, respectively. If the modulated signal r(t) is made 
periodic with period T 0 , then the complex envelope can be represented by the 
Fourier series 



f(t) = 




(106) 



where 



1 TJ2 

b k = — f r«) it , 

T„J „ 



■ o v - T /o 
o " 



(107) 



and f Q = 1/ T c . Equation (106) assumes a complex envelope bandlimited to Kf 0 . 
From Equation (106), the transform of the complex envelope is given by 

£ 

h/)= X ‘k *(/- • ( 10S > 

k = -£ 
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If the complex envelope is sampled with sampling period T s = T 0 /(2K + 1), which 
is basically the Nyquist rate, then, from Eqs. (106) and (107), 

[ 

f(lT s )= ^ b k e j2 * kl t L , (109) 

k = -/ 

and 



b k = - ^ r(lTs) e ~ j ' 2 * kl / L , (110) 

1 l = -f 



where T s = T 0 /L is the sampling period, and A = 2A'+1 is the total (odd) 
number of time and frequency samples. From Eqs. (103) and (108), the signal 
transform is given by 



ft 

R(f) = 0.5 £ { b k 6(f-f e - kfj + 4* <(-/- / c - i/J } , (111) 

k = -J 

from which the sampled signal can be derived as 

r(/T s ) = Re | h k <J 2jkl l L <J 2j fc lT s 

U = -/ 

(2) Subroutine Description. Figure 8 lists the pseudo-code of the 
subroutine ccq. The main input data for the subroutine are: 



( 112 ) 
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• Pulse length Tp , and pulse period (or pulse repetition period) T 0 . 

• Swept bandwidth and pulse amplitude A. 

• Maximum order of frequency component to be processed, A' ma x- The number of 
frequency components to be processed is (2 A’max + 1). 

• Sampling rate multiplier Sk, a term to be multiplied by the Nyquist rate in 
order to obtain the actual sampling rate. 

• Carrier frequency f c . 

• A flag test used to indicate option for plotting. 



Subroutine Ccq 

Compute data for reference signal: 

total number of samples L for the complex envelope /Eq. (114)/ 
sample period TS for the complex envelope /Eq. (113)/ 
phase— deviation constant DP /Eq. (105)/ 

Generate samples of ref . signal complex envelope /Eq. (102)/ 
Compute complex Fourier coefficients /Eq. (110)/ 

IF (flag test is TRUE) THEN 

Compute and plot transmitted signal /Eq. (116)/ 

END IF 

End Subroutine Ccq 



Fig. 8. Simplified Pseudo-Code for the Signal Generator Subroutine. 
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The total number of samples for the complex envelope is given by 



L=T 0 /T S , (113) 

where T 0 and T s are the pulse repetition period and complex envelope sampling 
period, respectively. The sampling frequency f s = 1/T S is the product of the 
sampling rate multiplier Sr by the Nyquist rate, twice the complex envelope 
bandwidth. In terms of the input parameters listed earlier, Eq. (113) can be 
rewritten as 



/ s 

/ A V 

L z I 0 Sr 2 ( \Sv\ + 2/Tp ) , (114) 

' v/ ' 

z M( envelope) 

which is the expression used in the subroutine. If the number of samples is less than 
the number of frequency components to be processed, (2 K m ax + 1), L is increased 
accordingly. If the computed value of L is even, the result is increased by one. The 
phase -deviation constant Df and the sampling period Ts are computed according to 
Eqs. (105) and (113), respectively. 

The subroutine proceeds to compute the samples of the complex 
envelope of the reference signal — using Eqs. (102), (104) and (105) — at the time 
instants t = / Ts . The Fourier coefficients are computed next, using Eq. (110). 

The signal to be transmitted is defined by multiplying the Fourier 
coefficient series { 6 k } — computed for the reference signal — by a window of length 



56 



(2 A max + 1)- The simple truncation, which corresponds to a rectangular window, is 
avoided because of the ringing it would cause on the waveform. Instead, the 
Hamming window 

w(n) = 0.54 + 0.46 cos( * n/K max ) , -A' max < n < A' max , (115) 

is used. Therefore, the transmitted signal and its complex envelope are given by 

{ finax \ 

£ U(k) b k J 2M ! L J 2 U‘TA , -A < / < A , (116) 

jr — _ If f 

* ~ A m ax 

and 

An ax 

f(lT s )= £ w{k)b k e j2tkl/L ,-K<l<K, (117) 

k = 

The transform of the transmitted signal is given, from Eq. (116), by 

^max 

XU) = 0.5 £ «<*) ( k* <(/- i - y„) + t* S(-J- i - t/o) } ■ (116) 

k = "fiaax 

While the sampling rate is high enough for the complex envelope 
representation, the signal r(<) usually requires a higher rate. In order to recover the 
signal from the Fourier coefficients, the series {6 k } must be padded with zeros, so 
that the modulated signal samples occur more than twice per carrier period. Both 



the subroutines ccq and calyce , to be described next, use six samples per carrier 
period for plotting purposes, for a total of 6f c T 0 samples on the interval 
(-T 0 /2 , T 0 / 2). This zero padding is accomplished by using the value L c = 6f c T 0 
instead of L in Eq. (116). 

6. Subroutine Calyce 

The subroutine calyce computes the output electrical signal, both in 
complex envelope y(t) and in band limited y(t ) forms. 

(1) Sampled Output Signal. The electrical signal at the 

output of the element ( m,n ) in the receiver array is given by Eq. (40c). 
Substitution of the expression for the input signal transform, Eq. (118), into Eq. 
(40c) yields 



,(IT S ) = 




max 



w{k) b , 



m + *4.0 e~ j2 'l ITs \ , (119) 



k — “In 



after performing the integration indicated in Eq. (40c), and where f 0 T s = 1 /L . 
Prom the definition of the complex envelope, Eq. (101), the sampled complex 
envelope of the output signal is given, by inspection of Eq. (119), by 



"max 

V, m (IT,)= £ M.k)b k HU c + kf^)e>' 2 M l L 



( 120 ) 



k = -i' 



max 



Equations (119) and (120) are both implemented in the subroutine calyce The 
observations made about the recovery of the signal from its Fourier coefficients — 
regarding the zero padding of Fourier coefficients — are also valid for y(t). In order 
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to have a correct reconstruction of y(t), the constant L in Eq. (119) is substituted 
by L c = 6^T 0 , corresponding to six samples per carrier period. 

(2) Subroutine Description. Figure 9 lists the pseudo-code of the 
subroutine calyce. The main input data for this subroutine, besides those listed for 
the subroutine ccq , are a three-dimensional array with the values of the overall 
transfer function for each frequency (J c + kj Q ) and position (m,n) in the receiver 
array, and the Fourier coefficients b k computed in the subroutine ccq. 

The computation of the total number of samples L and sampling 
period T s are as described in the previous section, in connection to the generation of 
the transmitted signal. The complex envelope samples y mn (lT s ) are computed 
according to Eq. (120). The signal samples y mn (lT s ) are computed using Eq. (119), 
with zero padding. 



Subroutine Calyce 
Compute : 

total number of samples L for the complex envelope 
sample period TS for the complex envelope 
Compute complex envelope samples /Eq. (120)/ 

Compute and plot signal /Eq. (119)/ 

End Subroutine Calyce 



Fig. 9. Simplified Pseudo-Code for the Subroutine Calyce. 
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2. Computing the Overall Transfer Function 

The computation of the overall transfer function is performed by a group of 
subprograms subordinated to the subroutine chfmn , called by the main program. 
The relationships between the subroutines used for this implementation are shown 
in Fig. 10. The program was implemented in FORTRAN, with most of the 
variables shared through common blocks. 

The pseudo-code for the subroutine chfmn is shown in Fig. 11. This 
subroutine implements Eq. (100), with the path of integration T truncated to the 
intervals (-jlower , j 0) plus (0 , upper). 

3. Integrand Implementation 

With reference to Fig. 10, the subprograms related to the computation of 
the integrand are reinteg , iminteg, integr, refl , hwvgky and dbsjO ’. They are all 
implemented as FORTRAN external functions. 

The functions reinteg and iminteg just return the real and imaginary parts 
of the integrand, respectively. They are called by the integration subroutine, whose 
first argument is the name of the function to be integrated. 

a. Function Integr 

The function integr computes the (complex) integrand in Eq.(100). 
The pseudo-code for this function is shown in Fig. 12. The argument k y 2 is 
interpreted as imaginary when it has a negative value, in all subprograms. The IF 
statements marked in the list are related to the medium transfer function. New 



■IMSL, Inc.. DBSJO , SFUNj Library. Bessel Functions of In teger Order , 1984. 
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Fig. 10. Structure of the Subroutine Chfrnn. 



G1 



Subroutine Chfmn 

FOR ( each frequency /= fc. + k X ) DO / ~^max ^ k ^ ^max/ 

EXECUTE Height /compute transmitter array weights (/)/ 

FOR C each position y n in receiver ) DO 

EXECUTE Break /compute breakpoints for integration/ 
Compute upper and lower limits of integration 
FOR ( each position x m in receiver ) DO 

EXECUTE QintegCReal part) /compute Re{ H(J, r mn ) } 
and its error Re{ EH(f,T im ) }/ 

EXECUTE Qinteg(Imag part) /compute lm{ //(/,r mn ) } 
and its error lm{ EH(f, r mn ) )/ 

Compute magnitude and phase of H and EH 
END FOR each x_ 

ill 

END FOR each y n 

Plot magnitude and phase of //versus position (x, p) 

END FOR each frequency 
End Subroutine Chfmn 



Fig. 11. Pseudo-Code for the Subroutine Chfmn. 
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Function IntegrC k y2 ) 

/this function is executed for constant frequency, receiver 
position ( m,n) and k y2 / 

Compute k r /Eq. (92)/ 

IF 1 (layered waveguide) THEN 

R 2l = Refl ( surf ace data) and /? 23 = Refl (bottom data) 

END IF 

/compute product ^-k 2 , for each y q , as an array Ftraniq )/ 

FOR ( each position y q in transmitter array ) DO 

IF 1 (layered waveguide) Ftran(q) = Hwvgky /k y2 •'K/, £ y2 . y Q \ y n )/ 
END FOR each y q 

/compute integrand according to Eq. (100)/ 

FOR (each position in transmitter array) DO 

Compute Bessel function J 0 (k r R 0 ) / R 0 is the horizontal 

distance from transmitter element to receiver position/ 
FOR ( each position y q in transmitter array ) DO 

Acumulate product c pq (/) •Ftran(q') • J 0 (k r R 0 ) / S S / 
END FOR each y q 
END FOR each r p 

Integr = (acumulated product)/27r 
End Function Integr 

■Other IF statements may be necessary for different medium transfer 
functions 



Fig. 12. Pseudo-Code for the Function Integr. 
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transfer functions are implemented by adding new IF statements. The functions 
integr and break to be described later, are the only ones that need to be modified 
when incorporating other medium transfer functions to the program. 

b. Function Refl 

This function computes the complex reflection coefficients. Its 
pseudo-code is shown in Fig. 13. 

The expressions for both reflection coefficients, Eqs. (73) and (76), can 
be written as 




i. . i. 

i’p *y2 “ Pi kyp 
f’p k y2 + Pi k yp 



(121) 



Function Refl(c p , a>, p p , p 2 , k. y2 ) 

/computes reflection coefficient /?2p/ 

Compute k. p /Eq. (92)/ 

IF (the speeds of sound are different) THEN 

Compute reflection coefficient by Eq. (l2l) 
ELSE 

Compute reflection coefficient by Eq. ( 12 2 ) 
END IF 

End Function Refl 



Fig. 13. Pseudo-Code for the Function Refl. 
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where k yp is given by Eq. (92). If the speed of sound is the same in both media, 
then Eq. (121) can be simplified to 




Pp ' Pi 
P p + Pi 



(122) 



c. Function Hwvghy 

The function hwvghy returns the complex product k y i #(/* y2>y 0 ;y n ) > 
the medium dependent factor in the integral expression for the overall transfer 
function, Eq. (100). This function is specific to the layered waveguide and 
implements Eq. (81), multiplied by k y2 . It is called by the function integr, as 
described earlier (see Figs. 10 and 12). Its implementation is trivial. Its arguments 
are the frequency, k y 2 and the y coordinates of the transmitter and receiver array 
elements, y and y n , respectively. The reflection coefficients and R 23 , and the 
depth D are available through common blocks. 
cL Examples 

In this section we examine the behavior of the integrand in the 
expression for the overall transfer function, Eq. (100). In all plots in which k y2 is 
the independent variable, negative values of k y . 2 are imaginary and positive values 
are real. For simplicity, both negative imaginary and positive real values of k y2 are 
plotted on the same axis. The symbol appearing on the plots is used by the 
plotting routine to mark different plots on the same graph, but it has no special 
meaning in the present context. 

(1) The Range Factor J 0 (k r R o ). Figures 14 and 15 show the plots of 
J 0 (k r R 0 ) vs Ay,, which is related to k r through Eq. (91). The frequency is / = 
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O.OG 




Fig. 14. Range Factor J 0 (k r R 0 ) for /= 500 Hz, Cj = 1,-500 m/s, and R 0 = 100 m. 
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Fig. 14. Range Factor J 0 (k r R 0 ) for /= 500 Hz, o, = 1,500 m/s, and 

R 0 = 100 m (cont.). 
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Fig. 15. flange Factor J 0 (k r R 0 ) for /= 500 Hz, = 1,500 m/s, and R Q — 1,000 m. 



68 



i-i 



OC 



0.8 



0.6 - 



0 . 4 - 



0.2 



0 - 



- 0.2 



MS 



^v«v«rr| 



0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 

KY2 (1/METER) 



Fig. 15. Range Factor J 0 {k T R 0 ) for /= 500 Hz, c, = 1500 m/s, and 

R 0 = 1000 m (cont.). 
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= 500 Hz and the sound speed is c^ — 1,500 m/s, which results in a wavenumber 
h ^ = 2.09 m ^ . In Fig. 14, the distance is R 0 = 100 m and, in Fig. 15, it is 1,000 m. 
In those figures the range of k y2 is (-y'0.25 , j 0 ) plus ( 0 , k^ ) . An important 

characteristic of the factor J o (k r R 0 ) is that it oscillates slower near k y2 z 0 . 
Another characteristic behavior is that it oscillates faster as the range increases. If 
the other terms in the integrand also have "slow" oscillations, then the greater 
contribution for the result of the integration comes from the region k y2 z 0 , as a 
result of the smoothing nature of the integration operation. Physically, the plane 
waves that travel near the horizontal, that is, those with k y2 ~ 0 , contribute more 
to the total field than those that interact more with the boundaries. In order to 
compare the "frequency of oscillations" of the Bessel function and of the medium 
transfer function, recall that, for k r R rj > 2 tt, the zeros of J 0 (k r R 0 ) occur 
approximately at intervals of 

6k=—. (123) 

R o 

(2) The Medium Factor. Figure 16 shows the plot of the term k y2 * 
for the layered waveguide whose characteristics are shown in Table 1. The depth D 
of the ocean is 100 m and the frequency is 500 Hz. The position of the point source 
is (jj, y s , Zs) = (0, 30, 0) meters and the point receiver sensor is at (.tv ,y r , £ r ) = 
= (1000 ,50 , 0) meters, corresponding to the range factor shown in Fig. 15. 

This waveguide is of the almost -pressure-release-surface/fast - 
bottom type. The peaks are evident in the figures. From the discussion in Section 
III -A -lb, there are 33 peaks in the range 0 < k y2 < 1.04 m ^ , separated by 
Ak y2 z ir/100. 
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Fig. 16. Waveguide Medium Transfer Function Times k y2 . (a) Real Part. 
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Fig. 16. Waveguide Medium Transfer Function Times k y2 . (a) Real Part(cont.). 
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Fig. 16. Waveguide Medium Transfer Function Times k yT (b) Imaginary Part. 
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Fig. 1C. Waveguide Medium Transfer Function Times fc y2 . (b) Imaginary 

Part(cont.). 
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TABLE 1. WAVEGUIDE ACOUSTIC CHARACTERISTICS 



Medium 


Speed of Sound 


Density 




(m s" 1 ) 


(Kg m 


I 


343.0 


1.21 


II 


1500.0 


1026.0 


III 


1730.0 


2070.0 



A criterion to set an upper limit of integration in Eq. (100) — the 
upper discussed in Section III-B-2 — could be based on the comparison between the 
"periodicities" of the medium and range factors in the integrand. If the factors in 
the expression for Eq. (81), are multiplied out, the faster oscillatory factor — the 
one with the greater exponent — has a "period" of tt/D or more, with respect to A\, 2 . 
that is, 




j 



D 



(124) 



From Eq. (91), the increments Sk y2 and Sk r are approximately related by 



\6k r \z 




(125) 



75 



when k y2 is real. Equations (123) and (125) can be used to compute the interval 
between zeros of the range factor, when it is plotted versus A; y2 . If we require, for 
example, that the integral be truncated when the Bessel function passes through 
zero jVtimes during one "period" of the medium factor, then, substituting Eq. (123) 
and Sk y2 = w / N D into Eq. (125), yields 

k 2 

upper = — ■ ■ ■■ — . ( 126) 

/ 1 + (Rol A r £) 2 

The upper limit computed according to Eq. (126) can smaller than the last peak 
r~r> j 

position, given by v k 2 - k 2 [see discussion following Eq. (99)]. In this case, the 
limit could be increased, or those peaks outside the limit could be simply discarded. 
In this implementation, the upper limit is set to ta- 

The lower limit of integration, in the region where £ y2 is 
imaginary negative and | k y2 • T | decreases exponentially as k y2 ~ -jx , should be 
set to minimize the error of integration for the slower decaying factor in k y2 'ty — the 
one with the exponent of smallest absolute value. The worst case is when there is 
one constant term, which happens when the source and the receiver are at the same 
depth, for example. Note that the term k^-ty is pure imaginary when k y2 is 
negative imaginary, and decreases towards zero as e = e *'^l^y 2 l. The lower 

limit is set to -j 0.125 k^ in this implementation. 

(3) The Integrand. Figure 17 shows the plot of the integrand for the 
waveguide discussed above, at / = 500 Hz, for a range Rq = 1,000 m. For & y2 
greater than about 1.5, corresponding to Nz 10 in Eq. (126), the integrand is fairly 
symmetrica] around zero, suggesting a low value for the integral in that region. 
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Fig. 17. Waveguide Integrand, Range Rq = 1,000 m. (a) Real Part. 
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Fig. 17. Waveguide Integrand, Range R 0 = 1,000 m. (a) Real Part(cont-). 
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Fig. 17. Waveguide Integrand, Range R 0 = 1,000 m. (b) Imaginary Part. 
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Fig. 17. Waveguide Integrand, Range R 0 — 1,000 m. (b) Imaginary Part(cont-). 
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4. Other Subroutines 
a. Subroutine Break 

This subroutine divides the interval of integration in order to decrease 
the error of the integration subroutine. The criteria used in this division axe the 
behavior of the medium transfer function and the relationship between the 
"periodicities" of that function and of the Bessel function in the integrand. For 
each different medium transfer function, a different interval subdivision strategy 
must be implemented. The upper and lower limits of integration are independent of 
the type of medium, in this implementation. Nevertheless, they could be made 
medium dependent, in which case their computation should be included in the 
subroutine break. Figure IS lists the pseudo-code for this subroutine. 

In the case of the layered waveguide function, the breakpoints are set 
equal to the locations of the peaks, if any, as discussed in Section III -A -lb, 
regarding the almost- pressure-release-surface/fast -bottom. The interval of values 
of k y2 outside the range of the peaks is divided evenly into subintervals of size 
bn / D , corresponding to five "periods," as given by Eq. (124). 

(1) Subroutine Description. The input parameters for the subroutine 
break are the frequency / the speed of sound c 2 — computed in the subroutine chjmn 
and taking into account the gradient of the speed of sound, for the case of other 
media. Other variables shared through common blocks are the type of medium, 
OCNTF, and the waveguide characteristics. The outputs are the number of 
breakpoints and a vector with the values of k y2 at the breakpoints. 

Initially, the subroutine divides the interval of integration evenly. 
Then, it proceeds to check for the presence of peaks in the medium transfer 
function, as discussed in Section III -A -lb. The peaks are the zeros of Eq. (99), 
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Subroutine Break( /, , breakpoints , number_brkpts) 

IF 1 (layered waveguide) THEN 

/divide range evenly (-0.125 k 2 , 0) U (0, kj/ 
number Jbrkpts = (1 .125 A£ y2 ) “ 1 

FOR in = 1 to numberjbrkpts ) DO 

breakpointsin) - (—0.125 hj + n A k y2 

IF ( breakpointsin ) < 0) index = n /index keeps track of 
last nonpositive element of vector breakpoints/ 

END FOR n 

/'check for presence of peaks/ 

IF (peaks do exist) THEN /a> rf2 , see Eq. (99)/ 

/the peaks are located at k y2 = 9 n -rD, see Eq. (99)/ 
numbe r p eaks - \ia/ir)- 0 . 5"j / discussion after Eq. (99)/ 

FOR in = inde. r+1 to indei+l+ numbe r p eaks ) DO 
Compute n_th_zero of the function Denwgd 
breakpoint sin) = n_th_zero -f Z) 

END FOR K 

Update numberjbrkpts to reflect total number of breakpoints 
END IF 
END IF 

End Subroutine Break 

'Other IF statements may be necessary for different medium transfer 
functions 



Fig. 18. Pseudo-Code for the Subroutine Break 
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implemented as the FORTRAN external function denwgd. The zeros are found by 
the subroutine dzbren 1 . The (positive) breakpoints originally computed are changed 
to the values of the positions of the peaks. This implementation assumes a lower 
limit of integration of -j 0.125 k^ and an upper limit equal to k 2 . 

b. Subroutine Weight 2 

The subroutine weight computes the transmitter array complex weights 
c p q (/). The input parameters are the frequency; the transmitter array data; the 
direction cosines u' (with respect to the X direction) and v' (with respect to the Y 
direction) of the transmitter array main lobe; the speed of sound at the center of the 
array — set to c 2 for the layered waveguide, but included in the present 
implementation to allow for other media — and the gradient of the speed of 
sound — set to zero in this implementation. The complex weights are given by [Ref. 

5:p.l33] 



c pq(f) ~ a pq 




(127) 



where 



CT 

CQ 

a 

II 

CT 

C? 


(128) 


$pq — -j- («' P d xt + v' q d yt ) , 


(129) 



•IMSL, Inc., DZBREN , MATH/ Library, Nonlinear Equations, 1985. 

2 This implementation is based on a program originally written by Dr. 
Lawrence J. Ziomek, at the Naval Postgraduate School. 



S3 



(130) 



A — ————— 

I 

A p and B q are the amplitude weights along the X and Y directions, respectively, A 
is the (depth -dependent) wavelength, d xt and d yt are the interelement spacing 
along the X and Y directions, respectively, and is the speed of sound at the 

depth y q . 

c. Subroutine Qinteg 

The subroutine qinteg is an interface for the subroutine dq2ags l , which 
actually performs the integration in Eq. (100). The input arguments are the name 
of the external FORTRAN function to be integrated, the limits of integration, the 
breakpoints, the number of breakpoints, and the required error in the result. The 
output arguments are the result of the integration and its estimated error. 

Initially, each interval of integration, defined by the breakpoints, is 
subdivided in order to improve the convergence of the subroutine dq2ags, which is 
called by qinteg to perform each intermediate integration. If the error of the 
integration of a subinterval, as computed by dq2ags, exceeds the requirements, that 
subinterval is further subdivided. Finally, the results and errors of integration 
provided by dqSags are accumulated. 

As indicated in the pseudo-code of the subroutine chfmn , Fig. 11, each 



•IMSL, Inc., DQ2AGS , MATH) Library, Integration and Differentiation , 

1985. 
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integration is performed twice, for the real and imaginary parts of the integrand. In 
terms of computing time, this is the most critical section of code. Not only does the 
subroutine qinteg calls dq2ags many times — at least ten times between breakpoints, 
in the present implementation — but, what is worse, dq2ags can execute the function 
to be integrated, ultimately integr, hundreds of times. The computation time 
increases with frequency and range: 

• The interval of integration increases as O(A^), directly proportional to 
frequency. 

• As the range increases, the range factor J 0 {k r R 0 ) oscillates faster, degrading the 
convergence of the integration subroutine, causing more executions of the 
function integr and more frequent interval subdivisions. 

The effect of range can be minimized if the limits of integration, lower and upper , 
decrease (in absolute value) with range, as in Eq. (126) for the upper limit. 

C. RESULTS 

The algorithms described in the previous section are tested by using two 
examples for which the results are easily interpreted. The unbounded homogeneous 
medium was presented in Section II-C-3a as a particular case of the layered 
waveguide, all three media with the same acoustic characteristics. The surface 
reflection problem was analyzed in Section II -C -3b. In both cases, the input to the 
system is a time-harmonic point source and the output electrical signal — 
numerically equal to the velocity potential — is given by Eqs. (85) and (88). From 
Eq. (40c), the time -independent term of the output signal is, for a time -harmonic 
input, equal to the overall transfer function, that is, 
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(40c) 



!U«)= )««r mn )e , ' 2r/t i/, 

-oo 

!W,(0 = J°KW "(A™) ^ 2,/t <V, (131a) 

-oo 

%nn ( ' 0 = H U, T mn) ^ 2 ^ ot ■ ( !3 lb ) 

Therefore, the magnitude of the overall transfer function is equal to the magnitude 
of the velocity potential. 

Results for the waveguide example used so far (see Table 1) are also presented. 

i u r\rr» /f ati nniin 
a. nuiuugcucuuo mcuiuiii 

Figure 19 presents the plot of the magnitude of the overall transfer function 
as a function of horizontal distance x from the source. The acoustic characteristics 
are those shown in Table 1 for medium II. The depth of the source is 30 m and of 
the receiver array, 50 m. The expected result is a straight, line on a log-log plot, 
corresponding to a velocity potential proportional to the inverse of distance. At 100 
m, the error of the computed value is 6%, and increases with distance up to 29% at 
3,100 m. The integration routine computed the integrand 3,360 times (total for 
real and imaginary parts) for 100 m of distance, a number that increased to 34,776 
at 3,100 m, indicating an increased difficulty in convergence of the integral. An 
analogous computation for the layered waveguide took 45,000 and 70,000 
computations of the integrand, for 100 m and 3,100 m, respectively. 
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Fig. 19. Overall Transfer Function for a Homogeneous Medium. 



2. Surface Reflection 

Figure 20 shows the overall transfer function as a function of depth, for a 
distance of 1,000 m between source and receiver. The acoustic characteristics of the 
two media are those shown in Table 1 for media I and II. The surface is, for all 
practical purposes, pressure release. The depth of the source is 30 m. For these 
data, the approximations for Eq. (89) are valid, indicating a sinusoidal variation 
with depth. The error increases from zero at the surface, up to a maximum error at 
30 m, the depth of the source — see discussion in Section III -B -3d, regarding the 
lower limit of integration. At 30 m, near a theoretical maximum, the error is 45% . 
At depths greater than 60 m, the error is negligible. For comparison purposes, the 
integrand was computed about 11,000 times for each value of the overall transfer 
function. 

3. Layered Waveguide: Waveform Prediction 

Figure 21 shows the transmitted signal. This signal was obtained by 
truncating the Fourier series of a CW pulse, as explained in Section III -B -la. The 
reference signal is a rectangular-envelope CW pulse of duration 20 ms and 
repetition period of 400 ms. The carrier frequency is = 500 Hz. Seventeen 
frequency components are used for the transmitted signal, that is, K max = 8. The 
resultant bell -shaped CW pulse has a duration of 100 ms. Ringing is barely visible, 
due to the low sidelobes of the Hamming window used for the truncation of the 
Fourier series. The extension of time in the plot of Fig. 21 is equal to one period. 
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Fig. 20. Overall Transfer Function for the Surface Reflection Problem. 
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Fig. 21 Transmitted Signal. 
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The acoustic characteristics of the media are the same as for the examples 
in Section III -B -3d and are listed in Table 1. The distance is R 0 — 1,000 m. The 
transmitter is located at (i^, y s , Zg) = (0, 30, 0) meters, and the receiver at 
(x r , Vri ^r) = (1000, 50, 0) meters. 

Figure 22 shows the received pulse when the transmitter is a single point 
source. The distortion due to the multipath nature of the medium is evident. 

Figure 23 shows the received pulse due to a 5*5 transmitter array. The 
interelement spacing is roughly A/2 in both the X and Y directions. Note that the 
pulse arrivals occur at the same instants as in Fig. 22, but the relative amplitudes 
are different. 

D. CONCLUSIONS 

The equations that describe a range -independent acoustic communication 
channel were derived by using linear systems theory, a basic engineering tool, as a 
framework. They incorporate aspects of signal processing and acoustic propagation. 
The main results from Chapter II are the expression for the overall transfer 
function, Eq. (51) or its alternative form, Eq. (60); and the procedure to obtain the 
transformed medium transfer function, the solution to the inhomogeneous 
Helmholtz wave equation. Eq. (90). These results, together with Eq. (40c). were 
used to implement a waveform prediction program, as described in Section III -B. 
We do not claim novelty of results. Equation (60) is a mere extension of a classical 
result, in which the contributions of the point sources that make up an array are 
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Fig. 22. Output Electrical Signal Due to a Single Point Source. 
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Fig. 23. Output Electrical Signal Due to a 5*5 Transmitter Array 
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simply added, a consequence of the linearity of the wave equation that describes 
small -amplitude acoustic phenomena. The medium transfer function for the layered 
waveguide was derived in order to enable us to implement a working computer 
program. 

The program described in Section III -B is based on a slightly modified form of 
Eq. (60), in which we use the wavenumber ^-component £ y2 as the independent 
variable. The only justification for this change of variable is the easier analysis of 
the medium transfer function. The program works, but is slow. The computation 
of the waveform shown in Fig. 23, for example, took about three minutes per 
frequency component, with an overhead time — independent of the number of 
frequency components — of about two minutes, in an IBM 3033U/4381 2-cpu 
network. Both the speed and accuracy can be improved by using an adaptive 
technique for the truncation of the interval of integration, that is, computation of 
the limits of integration as a function of the behavior of the medium and range 
factors in the integrand. 

Other implementations are possible, for example, the Fast -Field -Program 
described in [Ref. 4:pp. 90-92]. Equation (51). in the form of a two-dimensional 
Fourier transform, could be useful for implementation. 

Concerning the medium transfer function, analytical solutions to the 
range -independent wave equation are only available for a few speed -of -sound 
profiles. Approximate solutions, like those provided by the WKB method, are 
candidates for implementation. On the other hand, numerical methods may not be 
feasible for our purpose, waveform prediction, due to time constraints. 
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As mentioned in the beginning of Chapter III, the main thrust for the present 
implementation of the range-independent equations was to have a working program. 
The program is working. More than an improvement, it possibly needs a different 
implementation. Nevertheless, the basic important fact is that, independent of the 
implementation, the modular nature of the equations should be fully exploited. 
Solutions to specific problems can be implemented by adding new subprograms, 
with a minimum of code change. Further investigation is recommended. 
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